Color centers in NaCl by hybrid functionals 
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We present in this work the electronic structure and transition energies (both thermodynamic 
and optical) of CI vacancies in NaCl by hybrid density functionals. The underestimated transition 
energies by the semi-local functional inherited from the band gap problem are recovered by the PBEO 
hybrid functional through the non-local exact exchange, whose amount is adjusted to reproduce the 
experimental band gap. The hybrid functional also gives a better account of the lattice relaxation 
for the defect systems arising from the reduced self-interaction. On the other hand, the quantitative 
agreement with experimental vertical transition energy cannot be achieved with hybrid functionals 
due to the inaccurate descriptions of the ionization energies of the localized defect and the positions 
of the band edges. 
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I. INTRODUCTION 

The anion vacancy is a prototypical point defect in 
alkali halide crystals, which can be introduced by heat- 
ing them in an alkali metal vapor or by x-ray irradiation 
and by electron bombardment. Anion vacancies in alkali 
halide are color centers (F centers) because they are re- 
sponsible for the coloration of the otherwise transparent 
crystals^ The neutral anion monovacancy is the simplest 
form of the F center, with a single bound electron in the 
vacancy center. The electronic transition between the 
defect-induced levels can be trigged by the absorption of 
visible light. For example, the presence of the absorp- 
tion peak (F band) associated with the F center at 465 
nm renders the NaCl crystal yellow^— An eminent ap- 
plication of anion vacancies in alkali halides due to the 
coloration is the color center laser first demonstrated in 
1965 by Fritz and Menke^Z. 

Due to their exotic features and long history, F cen- 
ters in alkali halides have been studied extensively us- 
ing various techniques such as optical absorption^—, 
Raman spectroscopy^— , the Stark effect 1 ^!, and 
luminescence^—. It is now well established that the de- 
fect states induced by anion vacancies in alkali halides are 
localized deep levels because of the small dielectric con- 
stant, large band gap, and high effective mass of the host 
crystals. These localized defect states cannot be correctly 
described within the effective mass approximation. 21 
Early calculations employed either a linear combination 
of atomic orbital (LCAO) or the vacancy-centered vari- 
ational wavefunction method within the Hartree-Fock 
(HF) theory. The F center was treated as a hydro- 
genic system trapped in the Coulomb potential created 
by the ions as point charges^ - — Although these calcu- 
lations agree well with experimental data, it was later 
pointed out by Murrel and Tennyson, on basis of their 
LCAO self-consistent field (SCF) calculations within em- 
bedded cluster method, that the agreement could be for- 



tuitous, provided the results are sensitive to the selec- 
tion of model Coulomb potential surrounding the an- 
ion vacancy^ Thus it is desirable to use a "true" po- 
tential from ab initio calculations. Numerous LCAO 
SCF cluster calculations have been performed at the 
restricted open- shell (ROHF)2£i2 7 and unrestricted HF 
(UHF) level^r— The improvements, however, are lim- 
ited and not unanimous as a result of the lack of elec- 
tronic correlations and of the small cluster size limited 
by the computational capability. 

Recent advances in density functional theory (DFT) 
makes it the tool of choice for studying the electronic 
structure of point defects in insulators. One major defi- 
ciency of DFT is that the standard local and semi-local 
approximations for the exchange-correlation (XC) func- 
tional suffer from the spurious electron self-interaction 
and the lack of the derivative discontinuities of the XC 
potential with respect to the particle number—. The cal- 
culated band gap is therefore generally smaller than ex- 
perimental value, and the underestimation is more pro- 
nounced for wide band gap insulators. For NaCl, DFT 
calculations predict band gaps from 4.5 eV to 5.0 eV— ~— , 
much smaller than the experimental value (8.5 eV) 35 . 
The band-gap problem poses serious uncertainties to the 
positions of the defect electronic levels and makes the 
credibility of the calculation questionable. One of the 
simplest yet the crudest way to overcome the gap prob- 
lem is to apply a scissors operation to the conduction 
band and shift it rigidly to the experimental value rela- 
tive to the valence band maximum (VBM). However, this 
still leaves ambiguities to the position of defect levels. 
In particular, for deep levels which do not follow either 
the characteristics of the VBM or conduction band mini- 
mum (CBM), the scissors scheme is not well-justified and 
not satisfactory in first-principles calculations. Other 
physically intuitive approaches to the problem include 
self-interaction corrections (SICs) 36 , and the DFT+C/ 
method where the localized states with a strong Coulomb 
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repulsion, such as incomplete d or / shells in transition- 
metal oxides, are modeled by a Hubbard-like U term 37 
to open the band gap. 



An alternative, and a more general approach to the 
gap problem is to relieve the unphysical self-interaction 
and derivative discontinuity problems by incorporating 
the non-local exact exchange (or Fock exchange) into the 
DFT XC functional. One of the implementation of the 
exact exchange is the screened exchange (sX) method, 
which replaces the whole local-density approximation 
(LDA) exchange by a Thomas-Fermi screened Coulomb 
exchange potential 3 -^. In hybrid density functionals, a 
fraction of exact exchange is mixed with LDA or gen- 
eralized gradient approximation (GGA) exchange. The 
hybrid functionals generally not only improve the band 
gap, but also yield better results in bulk properties such 
as lattice parameters, bulk moduli and heats of formation 
for semiconductors and insulator a 39 ' 40 . A recent study 
comments that hybrid functionals are superior to SICs 
in reducing self-interaction errors since the unitary in- 
variance is preserved4i Successful applications of hybrid 
functionals in the defect properties of various oxides and 
semiconductors have been recently reported42r— 



In this work, we revisit the anion (CI) vacancy in NaCl 
in different charge states (-1, and +1) by both semi- 
local and hybrid functionals, aiming for a better under- 
standing of the F center in alkali halide and the perfor- 
mance of the hybrid functionals for localized defects. We 
note that besides the band-gap problem, the finite-size 
effect arises within the current modeling scheme for the 
vacancy. The common procedure for defect energetics 
calculation is to embed the defect into a supercell un- 
der periodic boundary condition (PBC). The advantage 
of using supercells instead of cluster methods is that the 
band structure of the host crystal is well-defined, as the 
cell is bulk-like. 49 However, tractable DFT calculations 
are usually constrained to about 1000 atoms, and the 
size of system is further limited for hybrid functionals in 
a plane-wave basis set. A single CI vacancy in a 1000- 
atom NaCl supercell corresponds to a vacancy concentra- 
tion of the order of 10 20 cm~ 3 , which is much higher than 
those found in experiment (10 15 to 10 19 cm -3 ). The pe- 
riodic images of the point defects in a high density thus 
give rise to unrealistic defect-defect interactions, making 
the defect energetics dependent on the size of the su- 
percells. The problem is even more serious for charged 
defects, as the neutralizing background slows the con- 
vergence of defect energies with respect to the supercell 
sizei££ Other sources of error for small supercells involve 
the elastic energy due to artificial relaxations of ions, 
and defect level dispersions introduced by defect-defect 
interactions. Corrections for the finite-size effect have 
been found indispensable in defect calculations for real- 
istic interpretations^—, and they will be discussed and 
applied to the present study of CI vacancies. 



TABLE I. Calculated lattice constant (ao), fundamental band 
gap at r (E g ), dielectric constant (too) and enthalpy of forma- 
tion (AHf) of rocksalt NaCl using various DFT functionals. 





a (A) 


E g (eV) 


Coo 


AH f (eV) 


GGA-PBE 


5.69 


5.00 


2.33 


-3.69 


HSE06 


5.65 


6.43 


2.13 


-3.85 


PBEO (a = 0.25) 


5.64 


7.19 


1.98 


-3.85 


PBEO (q = 0.40) 


5.62 


8.47 


1.86 


-3.93 


Expt. 


5.57^ 


8.5^ 


2.3^ 


-4.26^ 



a Reference 
b Reference 
c Reference 



II. CHOICE OF HYBRID FUNCTIONALS 

We first assess the performance of the hybrid function- 
als on the bulk properties of NaCl, and compare it to 
the GGA-PBE (Perdew-Burke-Ernzerhof) functional^. 
Both unscreened and screened hybrid functionals are em- 
ployed. In the unscreened PBEO functional, the exchange 
part of the XC energy E xc is constructed by mixing a 
fraction (a) of non-local exact exchange E x with PBE 
exchange E X BE , while the correlation energy is simply 
taken from PBE, £,f BE : 



(1 - a)E : 



PBE 



E, 



PBE 



(1) 



The amount of exact exchange a is a variable from to 
1, although conventionally a — 0.25 is used as suggested 
by perturbation theory^. In practice a is usually varied 
to meet the experimental gap value. In a plane-wave 
basis set, the evaluation of the exact (HF) exchange is 
a hog to the computational resources and tends to be 
rather slow because of its truly non-local nature. The 
calculation can be accelerated by truncating the slowly 
decaying long-range part of the exact exchange as in the 
HSE (Heyd-Scuseria-Ernzerhof) hybrid functional^: 

pHSE cnsiv.A , n „.i cnPBE,: 



E 



aE s J(v) + (1 
+^ PBE ' lr (M) - 



rf»BE 



(2) 



The screening parameter [i in Eq. ^ determines the sep- 
aration of the short-range (sr) and long-range (lr) parts: 

i = sr(r) + Ir(r) = 1 " + (3) 

r r r 

In one limit when /j, = 0, the long-range term is zero and 
HSE reduces to the unscreened PBEO functional. For 
/i — > oo, HSE is identical to GGA-PBE since the whole 
exact exchange is screened. Here we use the optimized 
fi = 0.207 A" 1 , following Ref. [H along with a = 0.25 
and refer to this functional as HSE06. 

In Table [J selected bulk properties of NaCl calculated 
using the GGA-PBE and hybrid functionals are summa- 
rized together with the experimental values. The calcu- 
lations are carried out in the projector augmented wave 
(PAW) framework with the VASP code^ir— . A semicore 



pseudopotential of Na is used, treating the 2p3s elec- 
trons as valence electrons. The kinetic cutoff energy for 
the plane- wave basis set is 500 eV. A T-centered 8x8x8 
Monkhorst-Pack k-point mesh^ 4 - is applied to the primi- 
tive cell containing one formula unit of NaCl. For HSE06 
and PBEO calculations, a down-sampled 4x4x4 mesh 
is used to evaluate the non-local exact exchange. The 
down-sampling for the non-local exchange reduces the 
computing time significantly. It is generally necessary 
for the PBEO to have a finer k-point mesh than for 
the screened HSE functional to reach convergence 
For the present case, the HF exchange using the PBEO 
changes by roughly 15 meV per atom from the down- 
sampled 4 x 4 x 4 to the full 8x8x8, while the energy 
is already converged within 10~ 2 meV with the HSE06 
functional. Nevertheless, the choice of the down-sampled 
k-point for the PBEO calculations is sufficient for the 
bulk properties. The lattice constant a is determined 
when the residual force is smaller than 5 me V/A. The 
high frequency macroscopic dielectric constant can 
be calculated within a GW scheme using the random- 
phase approximation (RPA) with local field effect in- 
cluded. Around 90 empty bands are used for calculating 
the dielectric constant. The dielectric constant will also 
be referred to later for the finite-size corrections. Finally, 
the formation energy AHf is obtained 

AHf = £ Na Cl( s ) - ^Na(s) - 2 E c\ 2 {g)- ( 4 ) 

In Eq. dU, £? Na ci( s ) is the total energy of bulk NaCl. 
£ Na ( s ) is the energy of bulk Na in a body centered cubic 
(bec) , which was optimized and calculated using the same 
k-point mesh and cutoff energy as the bulk NaCl. E C i 2 ( g ) 
refers to the energy of one gas phase CI2 molecule in a 
large tetragonal cell. 

One immediately observes that the hybrid functionals 
improve not only the direct band gap (Fis — > Ti) but also 
the lattice constant and heat of formation compared to 
the GGA-PBE calculation in Table U in agreement with 
earlier calculations . 39 i 40 Yet, it is found the band gaps 
are still underestimated for the hybrid functionals with 
the original fraction (0.25) of exact exchange, and the 
PBEO yields a much closer value to experiment than the 
HSE06. This implies that for wide gap insulators, as the 
electronic screening is quite weak, the unscreened exact 
exchange in PBEO is preferred. For defect calculations, it 
is customary to tune the fraction of the exact exchange so 
that the experimental band gap can be reproducedi 47 i 65 
By increasing the amount of the non-local exchange from 
0.25 to 0.40 within the PBEO, the band gap of NaCl re- 
covers nearly to the experimental value, and the lattice 
constant and heat of formation are also reproduced best 
among the chosen functionals. We note that hybrid func- 
tionals tend to underestimate the dielectric constant of 
NaCl, a trend also found for semiconductors and other 
insulators i££ An accurate description of the electronic di- 
electric constant with the hybrid functionals will require 
an explicit account of excitonic effects i 67 i 68 
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FIG. 1. (color online). The electron density (V° and V~) and 
hole density (V + ) isosurface of the a\ g state in the (100) plane 
calculated with the mPBEO functional. The lines are drawn 
in intervals of 0.01e/A 3 . The displacements of the nearest 
neighbor atoms around the CI vacancy after relaxation are 
illustrated by the arrows. The blue and green circles represent 
the Na and CI atoms, respectively. 

To this end, we face several functionals for the subse- 
quent calculations of the CI vacancy in NaCl. The PBEO 
(a = 0.40) (we will refer it to mPBEO hereafter) is appar- 
ently favored since it reproduces the experimental gap. 
However, as the choice of the fraction of the non-local 
exact exchange is empirical to some extent, its impact 
on the position of the deep defect level for wide gap insu- 
lators is still unknown. Meanwhile, the screened hybrid 
functional is of great interest as it shows considerable 
success in the prediction of defect energetics. Therefore 
it is plausible to also include the HSE06 functional with 
the original a, as well as the GGA-PBE for the defect 
calculations. 



III. ELECTRONIC STRUCTURE OF 
CHLORINE VACANCIES 

In this section we briefly sketch out the single-particle 
Kohn-Sham (KS) eigenvalues of the CI vacancy induced 
electronic levels. Supercells containing 64 atoms are em- 
ployed for the calculations. The Brillouin zone is sampled 
with a 2 x 2 x 2 k-point mesh, and a plane- wave basis set 
cutoff energy of 450 eV is used. Further, the k-point mesh 
for the non-local exact exchange is down-sampled to the 
T-point for HSE06, while a full 2x2x2 k-point mesh 
is necessary for well converged energies in PBEO calcu- 
lations. The convergence criterion for full relaxations is 
0.01 eV/A. 

The removal of one CI atom in a perfect NaCl crystal 
leaves a neutral vacancy V with one electron bound to 
the vacancy center. The localized nature of the unpaired 
electron can be clearly identified in the charge density 
isosurface shown in Fig. [TJ The Is characteristics of the 
wavefunction in the vacancy is contributed equally from 
the six neighboring Na atoms. The negligible displace- 
ment of the neighboring atoms around V° (see Table [TTJ) 
keeps the singly occupied a\ g level unshifted after relax- 
ation. In the +1 charge state V + , the a\ g state is unoc- 
cupied and the polaronic hole is trapped in the vacancy. 
The nearest neighbor Na atoms tend to relax away from 
the vacancy because of the positive electrostatic potential 
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TABLE II. Displacements (in the equilibrium bond length 
0.5ao) of the nearest-neighbor Na atoms around the CI va- 
cancy calculated using the 64-atom cell. The positive value 
represents an outward relaxation against the vacancy, and 
vice versa. The values obtained with the 216-atom cell using 
the GGA-PBE are shown in parentheses. 





GGA-PBE 


HSE06 


mPBEO 


v + 


+0.035 (+0.043) 


+0.035 


+0.036 


v° 


-0.000 (-0.000) 


-0.000 


-0.000 


v~ 


-0.039 (-0.040) 


-0.041 


-0.044 



band gap. We note that the absolute KS energies of the 
dig states [eKs(°ig) + e VBM] are roughly unaffected when 
going from semi-local to hybrid functionals. In hybrid 
functionals, the VBM is lowered by 0.9 eV (HSE06) and 
2.0 eV (mPBEO) with respect to the GGA-PBE as a re- 
sult of the reduced self-interaction for the CI 3p states. 
On the other hand, for rigid structures, the widening of 
the band gap in hybrid functionals tends to place the a\ g 
of V~ further away from the CBM compared to the semi- 
local functional, while the unoccupied a\ g state of V + is 
closer to the CBM when hybrid functionals are used. 



TABLE III. Energy levels (in eV) of the single-particle Kohn- 
Sham dig state of the CI vacancy in a cubic 64-atom cell 
referenced to the VBM. The energy is averaged over the BZ. 
For V°, the positions of the ai state in spin-up (occupied) 
and spin-down (unoccupied) channels are given. The absolute 
positions of the host band edges are also given. 







PBE 


HSE06 


mPBEO 




rigid 


relaxed 


rigid 


relaxed 


rigid 


relaxed 




3.74 


5.02 


4.90 


6.25 


6.81 


8.25 


V°(a\ B h 


4.09 


4.10 


4.83 


4.83 


5.30 


5.30 




4.85 


4.86 


6.13 


6.13 


8.12 


8.12 


5.04 


4.37 


6.07 


5.22 


6.87 


5.69 


6VBM 




-0.76 




-1.63 




-2.78 


£CBM 




4.24 




4.80 




5.69 



inside the vacancy. The outward relaxation delocalizes 
the a\ g (see Fig. [T]) , and shifts it to higher energy towards 
the CBM. In the -1 charge state V~ , the a\ g state be- 
comes doubly-occupied. Upon relaxation the six nearest 
Na atoms show inward displacement towards the vacancy 
(Fig. [TJ as a polaronic distortion. As a consequence, the 
two electrons are more localized inside the vacancy site, 
shifting the ay g state to lower energy. We note that the 
relaxations of the neighboring atoms around the anion 
vacancy follow the Oh symmetry for all charge states. No 
symmetry lowering (or Jahn- Teller distortion) is found 
since the defect level is either singly occupied for V°, or 
doubly occupied (unoccupied) for V~ (V + ). The reduced 
self-interaction in hybrid functionals also results in more 
pronounced atomic displacements for the charged defects 
as seen in Table [TTl 

Table Hill summarizes the KS energies of the CI vacan- 
cies in various charge states for both rigid and relaxed 
defect structures. The choice of k-point gives rise to 
a dispersion of the electronic level within the finite-size 
supercell scheme. The dispersion introduces a strong 
dependence on the supercell size of the energy level at 
the T-point e r . Thus we average the defect level en- 
ergy over the Brillouin zone e since the averaged level 
shows a much better convergence than e r £i The disper- 
sion also slightly pushes the host CBM to higher ener- 
gies. The finite-size effect will be discussed in detail in 
Sec. lIVBl As predicted by all functionals, the a\ g states 
in all charged states lie within the upper half of the host 



IV. THERMODYNAMIC TRANSITION 
ENERGIES AND FINITE-SIZE CORRECTIONS 

In general the single-particle energy level of the defect 
as calculated from the KS equation differs from the ex- 
perimentally observed transition energies^ A rigorous 
approach to the transition energies, as discussed in this 
section, relies on the total energy difference of the defect 
energetics in various charged states. 

A. Formalisms of formation energies and transition 
energies 

A central quantity for the defect energetics is the for- 
mation energy Ef for a defect D in charge state q 

Ef = Er> - En + BriifXi + g(e V BM + (5) 

where E-q and En are the total energy of the supercell 
with the defect D, and the host supercell without de- 
fects, respectively. n$ is the number of atoms removed 
from the supercell (jij — 1 for the CI monovacancy) or 
the number of impurities added (m < 0). /ij refers to 
the chemical potential of the associated defect particle 
reservoir, and is subject to equilibrium conditions. For 
the present study, under extreme Cl-rich (or equivalently 
Na-poor) conditions 

1 

Mci = 2 cl 2(g) - ( 6 ) 

This places an upper limit on fj,ci- The lower bound can 
be deduced from the following relation: 

MNa + MCI = ^NaCl(s)- (?) 

Therefore under Cl-poor (or Na-rich) conditions, which 
facilitate the formation of CI vacancies, 

MCI > ^NaCl(s) - MNa(s) = A iff + --Eci 2 (g)i (8) 

and this sets the lower limit of /ici- 

The remaining term £vbm + £ f m the formation en- 
ergy [Eq. ([5)] represents the chemical potential, or Fermi 
energy of the electrons in charged defects. The Fermi 
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energy €f is varied within the band gap referenced to 
the energy of the host VBM e V BM (0 < ep < E g ). Here 
£vbm is evaluated as the energy difference between a per- 
fect host supercell and the same host supercell with one 
electron removed from the VBM: 



e V BM = £°(n)-£+(n-l), 



(9) 



where n is the number of electron in the host supercell. 
In principle, one needs a sufficiently large supercell with 
n — > oo corresponding to the dilute limit. In practice, a 
fractional charge q can be used along with a small super- 
cell to obtain the £vbm 



£vbm 



lim -2 — 

q->0 q 



(10) 



In the present case, the cvbm converges well within a 
64-atom cell and a fraction charge of 0.001 e. 

For charged defects, it is evident from Eq. ((SJ) that 
the formation energy is dependent on the chemical po- 
tential of the exchanged electron. The thermodynamic 
transition energy e(q/q') is defined as the Fermi energy 
at which the charge state q and q' of the defect sys- 
tem can be transformed spontaneously from one to the 
other. Therefore at the transition energy e(q/q') these 
two charge states have the same formation energy. This 
gives the following form of the transition energy 



<q/q') 



E u (q)-E D (q') 



£vbm- 



(11) 



B. Finite-size corrections 

Before proceeding to the results, we shall discuss the 
correction methods for the finite-size effect, whose causes 
have been already identified in Introduction. For charged 
defects, the simplest correction is to align the electro- 
static potential in the defect supercell to that of the host 
supercell. This is usually done by inspecting the potential 
difference AV between the core potentials of the atoms 
far from the defect center and that of the bulk cell, and 
the energy correction term is essentially AE = qAV. 
This correction is rationalized by the fact that in periodic 
supercell calculations the zero of the electrostatic poten- 
tial is chosen arbitrarily for each calculation, and the 
charged defect gives rise to a constant shift in the poten- 
tial so that the bulk VBM cannot be applied directly to 
the defect supercell. However, due to the small size of the 
supercell used even the atoms farthest from the charged 
defect center are not bulk-like, making such correction 
scheme inaccurate. Recent study reveals that the po- 
tential alignment resembles the Makov-Payne scheme^ 
whereas the latter targets the correction of the unphysical 
defect-defect interactions. Indeed, Komsa and Rantala 
found that AV has the form (eL) _1 , where L is the lat- 
tice constant of a cubic supercell^ This is analogous to 
the Makov-Payne scheme to first order. 



The popular Makov-Payne scheme for a charged de- 
fect in a cubic supercell in the dilute limit (L — > oo) is 
expanded as 



Ef{L) — Ef(L 



2eL 



2irqQ 
3eL 3 



+ 0{L~% (12) 



where «Md is the Madelung constant dependent on the 
lattice type, and q and Q the monopole and quadrupole 
moment of the defect charge, respectively^. The first 
order term, also called the Madelung energy, is thus the 
correction to the monopole-monopole interaction arising 
from the periodic image. The L _1 behavior of the arti- 
ficial electrostatic interaction vanishes slowly, and this is 
usually the leading source of error. The higher order cor- 
rections have much smaller effects on the formation en- 
ergy for ionic crystals, and it is usually accurate enough 
to include the corrections up to the quadrupole term. We 
note that in the Makov-Payne scheme the defect states 
are assumed to be localized, which is the case for the 
CI vacancies in NaCl. For delocalized levels higher or- 
der corrections might become necessary. Although the 
Makov-Payne expansion is sound and accurate, it has 
been found that direct corrections using the Madelung 
energy and multipole interactions are prone to overshoot 
the formation energy, in particular for small supercellsZS. 
A more reliable approach is to employ a scaling method 
by performing a series of calculations using supercells of 
different sizes with the same symmetry^ The corrected 
formation energy Ef{L — > oo) then can be extrapolated 
to the dilute limit by fitting the calculated formation en- 
ergies within finite-size cells to 



Ef(L) = Ef(L ->■ oo) + a\L 1 + a 3 L 



-3 



(13) 



where a n and Ej(L — ¥ oo) are fitting parameters. It 
is clear that this scaling law method requires at least 4 
supercells, and is rather computationally laborious. 

In a recent work Freysoldt et al. proposed a general 
correction scheme (we will refer to it as the FNV scheme 
hereafter) for finite-size effect based on a single calcula- 
tion of defect supercell without empirical parameters^: 



E f = Ef(L 



oo) + 4 att 



q\/b, 



(14) 



where £^ att is the macroscopically screened lattice energy 
of the defect charge qd with compensating background, 
and A q / b is an alignment term referenced to the bulk 
supercell to account for the microscopic screening. As 
the long-range ii^ att scales as L^ 1 and the short-range 
alignment term as L~ 3 , the FNV scheme can be seen 
as an extension to the Makov-Payne expansion. It also 
allows for an explicit expression for the third-order L~ 3 
energy term. 

Now we apply both Makov-Payne scaling and FNV 
schemes to the formation energies of the charged CI va- 
cancies (V + and V~) in NaCl. We refrain from including 
the potential alignment in these two schemes in order to 
avoid double-counting of the long-range L _1 term. In- 
deed Castleton et al. noticed that finite-size scaling with 
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FIG. 2. (color online). Demonstration of the correction 
schemes for the formation energies of the CI vacancy in the +1 
and —1 charge states (in Cl-rich limit) with respect to the re- 
ciprocal supercell lattice constant L~ x . The calculations were 
performed with the GGA-PBE functional without structural 
relaxations. 



potential alignment resulted in wide error barSf^l A sc- 
ries of simple cubic supercells containing 64, 216, 512 and 
1000 atoms is chosen in the present study. The exceed- 
ingly large 1000-atom supercell restricts the calculations 
to the GGA-PBE functional, although we show that the 
obtained trend is applicable to hybrid functionals as well. 
For the 64- and 216-atom cells, the Brillouin zone is sam- 
pled with a 2 x 2 x 2 k-point mesh. For larger supercells, 
we use two special k-points, i.e. T-point (0,0,0) and R- 
point (0.5,0.5,0.5) in reciprocal coordinates. Moreover, 
in the FNV scheme, the point charge qd consists of an 
exponential decaying term and a localized contribution 
modeled by a Gaussian 



q d {r) = qxN ie - r ^ + q(l - x)N p e 



(15) 



where jV 7 and Np are normalization constants, and x 
is the fraction of the relative amount of the exponen- 
tial decay. In practice, the resulting corrected energy 
is insensitive to the choice of the specific parameters in 

Eq. C5j,2a 

In Fig. [5J we demonstrate the effects of finite-size cor- 
rections to the formation energies of CI vacancies in +1 
{V + ) and -1 (V~) charge states. No relaxation is taken 
into account at this stage so as to exclude the finite- 
size effect of elastic energies. For V + , one first notices 
that the extrapolated formation energy from the Makov- 
Payne scaling law falls in line with that of the FNV 
scheme. The L~ l clearly dominates for the Makov-Payne 
fitted curve. The FNV scheme, on the other hand, shows 
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FIG. 3. (color online). Band structures (upper panel) of 
the CI vacancy in the —1 charge state (V~) in the unrelaxed 
64- and 216-atom cells using the GGA-PBE functional. The 
shaded areas show the bandwidth of the doubly-occupied ai g 
state induced by the negatively charged defect. The horizon- 
tal dashed line indicates the CBM of the perfect host without 
defects. In the bottom panel, the evolution of the bandwidth 
of the a\ g state (by closed circles •) as a function of the super- 
cell lattice constant is shown with an exponential fit curve for 
the unrelaxed structure. The band dispersion for the relaxed 
cell is also given by open circles o. 



a rapid convergence of the V + formation energy. We see 
that finite-size correction is indeed mandatory for an ac- 
curate description of formation energy of charged defects. 
For the smallest 64-atom cell, the uncorrected formation 
energy is underestimated by roughly 0.6 eV. Even for the 
1000-atom cell, the formation energy without correction 
is still 0.2 eV too low. 

Complexity arises when we move to the -1 charged CI 
vacancy. The uncorrected formation energies in Fig. |2] for 
V~ exhibit a zigzag evolution with respect to the increas- 
ing supercell size, making the Makov-Payne fit unreliable. 
Meanwhile, the FNV correction apparently yields a too 
high energy for small supercells, and the value does not 
appear to converge until we use the 512-atom cell. The 
source of such error is identified as the spurious disper- 
sion of the defect levels as a result of the overlap between 
the wavefunctions of the defect and its periodic images. 
In the dilute limit, this localized defect level should be 
strictly a flat band. However, as shown in Fig. [3J the 
CI vacancy induced a\ g level within the band gap shows 
a prominent dispersion for small supercells. For hybrid 
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functionals, the dispersion is less pronounced as the ex- 
act exchange favors a more localized electronic state in 
the vacancy. This artificial interaction tends to push the 
CBM of the host crystal to higher energies, or otherwise 
the CBM will become populated. The defect level dis- 
persion has a short-range characteristics, as is evidenced 
by the exponential fit of the ai g bandwidth with respect 
to the supercell lattice constant L for the 64-, 216- and 
512-atom cells (in Fig. [3]). This short-ranged effect is 
not included in either the Makov-Payne or FNV scheme, 
and thus one has to take it into account explicitly. We 
note that the dispersion correction is not necessary for 
V + since its a\ g level is unoccupied. Here the correction 
for the dispersion is considered by calculating the energy 
difference between the KS energy of the ai g level at the 
T-point e r and the ai g KS energy averaged over the sam- 
pled k-points in the Brillouin zone e. The FNV scheme 
based on the dispersion corrected formation energies is 
again able to yield converged results for small supercells, 
and the results are also comparable to the Makov-Payne 
method including the dispersion effect in the dilute limit. 

For the neutral CI vacancy V°, the situation becomes 
trouble-free since the electron is tightly bound to the va- 
cancy center with a strongly localized electron density 
distribution as seen in Fig. [TJ The formation energy 
barely varies for the supercells considered, and therefore 
there is no need for finite-size corrections for V . 

We have shown the finite-size corrections for the 
charged defect supercells in rigid geometries with atoms 
fixed at their bulk positions. However, the introduction 
of a CI vacancy inevitably changes the electrostatic po- 
tential of the local environment, resulting in atomic re- 
laxations around the vacancy. The supercell approach, 
in this aspect, will lead to another error because the su- 
percell employed in practice is usually not large enough 
for all local relaxations around the defect. This error can 
be partially alleviated by restricting the relaxations to 
the first two atomic shells around the defect, although 
it might underestimate the relaxation energy. Here we 
assess the finite-size effect on the elastic energy of the CI 
vacancy in various charge states based on full relaxations 
using the GGA-PBE functional. We do not discuss the 
relaxations of the outer shell atoms since their displace- 
ments are much smaller than the first shell Na atoms 
and they contribute little to the formation energy. For 
the neutral vacancy V , negligible inward relaxations of 
the six nearest neighbor Na atoms are found with super- 
cells containing up to 216 atoms (see Table |H|) . For V~, 
the net negative potential induced by the excess electron 
added to the vacancy gives rise to an inward displacement 
of the neighboring cations. As discussed in Sec. IIII1 this 
results in a more localized a\ g state, which consequently 
suppresses the dispersion of the a\ g with respect to the 
rigid structure (see Fig. [3]). Due to the finite-size of the 
supercell, Table |H] shows that the displacement obtained 
from the 64-atom cell is 0.05 A smaller than that from 
the 216-atom cell. The inability to fully relax in the 64- 
atom cell consequently gives a formation energy about 



0.2 eV higher than that of the larger supercells. For the 
positively charged vacancy V + , the six nearest neigh- 
bor Na atoms experience an outward displacement due 
to the positive potential in the vacancy center. In con- 
trast to V~, the finite-size effect on the elastic energy 
is not significant for V + , as seen in Table QT] since the 
atomic displacements using the 64- and 216-atom cells 
are of similar magnitude. 

With all these comprehensive finite-size effects in mind, 
we now summarize the correction scheme applied in the 
present work. We restrict the calculations of formation 
energies with hybrid functionals to the use of the 64-atom 
supercell. Thanks to the localized nature of the defect 
state and negligible relaxation for V°, no correction is 
necessary. For V + we apply the FNV to the 64-atom cell 
and refrain from any correction for the elastic energy. 
For V~ the dispersion correction is applied to the 64- 
atom cell, followed by the FNV scheme. Further, the 
relaxation energies are aligned with those obtained from 
the 216-atom cell, provided the latter already yields a 
converged elastic energy. In practice, due to the similar 
amount of atomic displacement (see Table |TT|, we use the 
PBE result as a reference for the hybrid functional, and 
subsequently lower the formation energies by 0.2 eV for 
the relaxed 64-atom supercells in the -1 charge state. 

C. Chlorine vacancy thermodynamic transition 
energies 

The calculated formation energies with finite-size cor- 
rections for the CI vacancy in NaCl are shown in Ta- 
ble IIV1 with the Fermi energy ep fixed at the VBM. 
For the neutral vacancy V , the hybrid functionals yield 
higher formation energies than the GGA-PBE, although 
the energy differences are small. The functional depen- 
dence of formation energies for the charged states is much 
more prominent. With the hybrid functionals, we obtain 
higher formation energies for V" and lower formation 
energies for V + . In particular, the negative formation 
energy of V + under Na-rich conditions suggests that the 
F + center could be predominating when the Fermi en- 
ergy is close to the VBM. 

To trace the source of the functional dependence of for- 
mation energy for charged vacancies, we may first rewrite 
the formation energy of V + in a rigid geometry with eF 
fixed at the VBM as 

E f (V+) = (E+ - Eft) + (£° - £° + M ci) + e VBM , (16) 

where the first term on the right-hand side is the electron 
ionization energy of V° (or equivalently the affinity en- 
ergy of V + ), and the second term simply the formation 
energy of V . By decomposing the formation energy into 
several contributions, it is clear that the discrepancies in 
Ef(V + ) stem mostly from the different positions of the 
VBM by various functionals. We note that the ionization 
energy of V° shows very small changes (within 0.05 eV) 
from semi-local to hybrid functionals, consistent with the 
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TABLE IV. Formation energies (in eV) of CI vacancies in various charge states calculated with the GGA-PBE, HSE06 and 
mPBEO functionals under Cl-rich and Cl-poor conditions. The Fermi energy is chosen at the VBM for charged defects. All 
values are corrected for finite-size effect. 







Cl-rich conditions 






Cl-poor conditions 




GGA-PBE 


HSE06 


mPBEO 


GGA-PBE 


HSE06 


mPBEO 


v u 














rigid 


4.44 


4.63 


4.71 


0.75 


0.78 


0.78 


relaxed 


4.44 


4.63 


4.71 


0.75 


0.78 


0.78 


V+ 














rigid 


1.47 


0.75 


-0.26 


-2.22 


-3.10 


-4.19 


relaxed 


0.63 


-0.10 


-1.10 


-3.05 


-3.95 


-5.03 


v- 














rigid 


9.71 


11.40 


12.91 


6.02 


7.55 


8.98 


relaxed 


9.55 


11.01 


12.50 


5.86 


7.16 


8.57 



similar absolute energy of the singly occupied a\ g state 
(see Table |ni|. Analogously, the formation energy of V~ 
can be rewritten as 



Ef(V) = (#d 



(££-£&+/zci)-evB M) (17) 



where the first term on the right-hand side is now the 
(negative) electron affinity energy of V° (or the ioniza- 
tion energy of V - ). In contrast to the ionization energy, 
it is found that hybrid functionals tend to yield a smaller 
affinity energy of V° than semi- local functionals, and that 
the difference reaches up to 0.9 eV. Along with the cvbm, 
they explain the variations in the formation energy ob- 
served in Table ITVl 

The atomic relaxation energy due to the polaronic elec- 
tron or hole can be further extracted from Table IIVI as 
the difference between the rigid and relaxed structures. 
It is not surprising that the relaxation energies given by 
various functionals are consistent provided the atomic 
displacements are similar with these functionals (see Ta- 
ble |TTJ) . The relaxation energy for V + is about 0.8 eV, 
while it ranges from 0.2 to 0.4 eV for V~. 

For a charged defect, the formation energy is a function 
of the Fermi energy as illustrated in Fig. @] The inter- 
sections of different charge states are the thermodynamic 
transition levels defined in Eq. (jlip . We see in Fig.@]that 
the transition levels e(+/0) and e(0/— ) are both within 
the band gap. Therefore, the mPBEO functional predicts 
that all charge states (-1,0 and +1) of CI vacancy could 
be thermodynamically stable when the Fermi energy is 
varied within the band gap. We note that although this 
is also qualitatively predicted by the uncorrected forma- 
tion energies, the Fermi energy window for the neutral 
V vacancy is much narrower. 

The density functional dependence of the thermody- 
namic transition levels is illustrated in Fig. EJa). Both 
neutral and +1 charge states are predicted to be stable 
since the e(+/0) levels are within the band gap for all 
functionals, with transition energies increasing from 3.81 
eV to 5.36 eV as the band gap widens from semi-local 
to hybrid functionals. On the other hand, the e(0/— ) 
is placed slightly above the CBM in the PBE, while its 
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FIG. 4. (color online). Calculated mPBEO formation energies 
of the CI vacancy with full relaxations as a function of the 
Fermi energy under the Cl-poor condition. The solid lines 
represent the formation energies corrected for the finite-size 
effect. The thermodynamic transition levels and the zero- 
phonon lines (ZPL) are indicated. The uncorrected values 
are given in dashed lines for reference. 
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FIG. 5. (color online), (a) The thermodynamic transition 
levels of the CI vacancy calculated with various functionals. 
The position of the VBM is aligned to energy zero, (b) The 
absolute values of the thermodynamic transition energies. 
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FIG. 6. (color online). Configuration coordinate diagram for 
the neutral and + charge state of CI vacancy in NaCl. The 
Fermi energy is located at the CBM. The optical processes 
involved are the absorption (a), emission (c), Stokes (b) and 
anti-Stokes (d) shifts, and zero-phonon line (ZPL). The coor- 
dinates at the lowest vibrational state of the ground state and 
excited state are denoted by Ro and R+, respectively. The 
zero-point energy is neglected in the diagram. 



position falls into the band gap in HSE06 and is fur- 
ther shifted downwards with respect to the CBM in the 
mPBEO calculations. Hence in contrast to the GGA- 
PBE, both the HSE06 and mPBEO imply that the V~ is 
stable. 

While the thermodynamic transition energy generally 
increases with respect to the VBM as the band gap en- 
larges, we see from Fig. [5jb) that the absolute position 
[e(q/q') + £vbm] of the e(+/0) level remains roughly un- 
affected from semi-local to hybrid functionals. This co- 
incides with the findings by Alkauskas et al. that the 
calculated energy levels of localized defect are generally 
not tied to the position of the CBMZL The e(0/-) levels 
are nevertheless more dispersed. 

V. OPTICAL PROPERTIES OF THE COLOR 
CENTER 

The experimentally available optical properties of the 
F (V°) and F' (V~) center in NaCl serve as a benchmark 
for the assessment of the performance of the functionals. 
The optical processes are clearly marked in the config- 
uration coordinate diagram in Fig. [6] according to the 
Franck-Condon principle. 74 In the Franck-Condon ap- 
proximation the electronic transition is assumed to oc- 
cur very fast compared with the motion of nuclei in the 
lattice. Therefore, the optical excitation spectrum ob- 
served in experiment does not involve the relaxation of 
the defect structure, in contrast to the thermodynamic 
transition. 

The optical absorption and emission can be described 
by vibronic (simultaneous vibrational and electronic) 
transitions, in which the lattice vibration mode is treated 
by a quantum harmonic oscillator. We first consider 
the excitation of an F center, which is a well-defined 
featured— ^ By absorption of a photon, the unpaired 



electron is transferred to an excited electronic state (V + 
state) and an excited vibrational state. The excitation 
of an electron into the CBM is equivalent to bringing 
an electron to a reservoir with a chemical potential of 
evbm + Eg. The optical absorption energy E a , as illus- 
trated in Fig. [BJ is thus given by 

E, = Ef°(+ ] e F = E g )-Ef°{0) 1 (18) 

where the first term on the right-hand side is the forma- 
tion energy of the unrelaxed +1 charge state (at coor- 
dinate Ro) with the Fermi energy at the CBM, and the 
second term the formation energy of the relaxed neutral 
charge state (at Ro). The excited state (F + center) sub- 
sequently relaxes to its zero-point vibration states. The 
energy gain due to the relaxation is the Stokes shift E$ 
between the vertical absorption and the zero-phonon line 
(ZPL). The ZPL is the transition energy from the lowest 
vibrational level (zero phonon mode) of the ground state 
to the lowest level of the excited state (at Ri), without 
energy transfer to lattice phonons. In terms of thermo- 
dynamic transition energy, it is easy to see from Fig. 0] 
that the ZPL can be expressed in terms of the difference 
between the band gap E g and e(+/0). In the present 
case, due to the identical formation energy between the 
rigid and relaxed V°, the Stokes shift reduces to the re- 
laxation energy of the V + from the rigid structure. In 
the vertical emission (luminescence) , the excited electron 
from the CBM recombines into the defect level and the 
emission energy E e is given by 

E e = E? 1 (+;e F = E g )-Ef 1 (0), (19) 

where the defect structure in the neutral state is kept 
fixed as that in the relaxed +1 charge state. Finite-size 
correction on this relaxation energy is taken into account 
according to Sec. UVBl Once the electron is in the ground 
state, it relaxes to the bottom of the state with the relax- 
ation energy Eas (*-e- the anti-Stokes shift between the 
ZPL and the vertical emission). The pronounced Stokes 
and anti-Stokes shifts (see Table fV|) are expected due to 
the large polaronic distortion. 

The calculated vertical absorption and emission ener- 
gies of the F center using the GGA-PBE and hybrid func- 
tionals are reported in Table [V] The zero-point energy 
is not included since it is usually comparable for both 
the ground state and excited state. It should be borne 
in mind that the total energy difference scheme (ASCF) 
based on the one-particle picture is not capable of de- 
scribing the electron-hole coupling in the optical absorp- 
tion. The presence of the exciton will usually introduce 
a uniform rcdshift of the absorption peak irrespective of 
the density functionals. 

We find that the absorption and emission energies are 
underestimated by the GGA-PBE, and this is most likely 
related to the small band gap. At the present stage it re- 
mains unclear whether mPBEO or HSE06 is more appro- 
priate for the absorption energy since the excitonic bind- 
ing energy of the F center is yet unknown. Meanwhile, 
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TABLE V. Calculated vertical absorption (E a ) and emission 
(E e ) energies, zero-phonon line (ZPL) and the Stokes (Es) 
and anti-Stokes (Eas) shifts of the F and F' centers in NaCl. 
All values in eV. 







£ a 


E c 


ZPL 




Eas 


F center 














GGA-PBE 




2.03 


0.65 


1.19 


0.84 


0.55 


HSE06 




2.56 


1.07 


1.70 


0.85 


0.63 


mPBEO 




3.50 


1.88 


2.66 


0.84 


0.78 


Expt. 




2.77^ 


0.98^ 








F' center 














GGA-PBE 




0.76 


-0.27 


-0.11 


0.87 


0.16 


HSE06 




0.91 


-0.34 


0.04 


0.86 


0.38 


mPBEO 




2.03 


0.27 


0.68 


1.35 


0.41 


Expt. 




2.43^ 










a Reference 
b Reference 
c Reference 


76. 
77. 



the available experimental ZPLs for several -F-aggregated 
centers (1.96 eV for R 2 bandZ 8 - and 1.48 eV for N band^) 
suggest the mPBEO might give a too large ZPL. Never- 
theless, the hybrid functionals yield more realistic optical 
transition energies than the GGA functional. 

Analogously, we extend the calculation to the optical 
process of the F' center, which is formed when an elec- 
tron is trapped at an F center by light absorption at 
low temperatures. - Instead of the sharp and bell-shaped 
curve of the F band, the F' center of NaCl gives rise 
to a broad F' absorption band, peaking at the longer 
wavelengths side of the F bandar— It is seen in Table IVl 
that all functionals now predict smaller absorption en- 
ergies with respect to the experimental F' band peak. 
This tendency is not changed even if the excitonic ef- 
fect is taken into account as the electron-hole interaction 
will further decrease the absorption energy. In accord 
with the F band absorption, the F' E a increases from 
the semi-local functional to the hybrid functionals as the 
calculated band gap widens. The E a values given by the 
GGA-PBE and HSE06 are well below the experimentally 
observed peak, whereas the mPBEO yields a value that 
is in better agreement with experiment. In addition, we 
find that the various functionals predict either a negative 
or a very small emission energy E c from the excited F' 
state to the ground state. A negative emission energy in 
Table \V\ suggests that the configuration coordinate of the 
intersection lies between the coordinates of the minimum 
of the ground state and excited state. It is conceivable 
that in such case the excited state can return back to 
the ground state through a non-radiative process by vi- 
brational relaxations, which leads to the luminescence 
quenching. The non-radiative path is also valid for a vi- 
bronic system with a small emission energy where the 
intersection is in the vicinity of the minimum of the ex- 
cited stated Experimentally, a radiative transition of an 
F' excited state is indeed absent^ 



VI. DISCUSSION 

We have shown that, while hybrid functionals have 
been reported to be adequate for defects in some 
semiconductors ^2r— the description of the localized an- 
ion vacancy in a wide gap insulator is less satisfactory 
by hybrid functionals when compared to the experimen- 
tal optical absorption spectra. In this section, we aim to 
identify the possible origins of the failure of hybrid func- 
tionals for the description of the color centers in NaCl. 

We start with the discussion of the optical absorption 
since it is well-defined experimentally. The absorption 
energy of an F center in Eq. (|18[) can be rewritten as 

E a = -e(+/0) Ro + E g 

= [£ R °(+) - £ R °(0)] + (6 V BM + E g ), (20) 

S v ' S v ' 

IP(V°) e CBM 

where e(+/0) R ° refers to the vertical transition energy 
occurring at the geometry for the neutral vacancy, and 
IP(V°) is the ionization energy of V°. Therefore the ab- 
sorption energy is solely dependent on the ionization of 
the neutral vacancy and the position of the band edge 
in a perfect supercell, and no structural relaxation is in- 
volved. As the ionization energy is not sensitive to the 
choice of the functional as discussed in Sec. IIVQ it be- 
comes obvious that the discrepancies in the absorption 
energy reported in Table [V] mostly stem from the varia- 
tions in the CBM energy. For instance, the GGA-PBE 
ccbm is 1-45 eV lower than the mPBEO value as a result 
of the well-known band gap problem associated with the 
local and semi-local DFT functionals. While the energy 
gap can be reproduced by mixing 40% exact exchange in 
mPBEO, it is yet not clear whether the positions of the 
band edges are accurate. 

In principle, the band gap problem can be overcome 
by quasiparticle (QP) self-energy calculations based on 
many-body perturbation theory. Here we follow the 
widely adopted GW approximation for the electronic self- 
energy 82 and calculate the QP corrections to the Kohn- 
Sham eigenvalues. The GW approximation can be un- 
derstood as the Hartree-Fock theory with a dynamically 
screened Coulomb interaction. The QP energies are cal- 
culated in a two-atom NaCl unit cell with a T centered 
4x4x4 k-point mesh and an energy cutoff of 200 eV for 
the response function, and a total of 256 bands. The dy- 
namic dielectric matrix is constructed with a frequency 
grid of 200 points^ We note that the QP gap of NaCl 
is closely related to the starting wavefunction and self- 
consistency. It is found that single shot GqWq correc- 
tion is too small when it is applied to the GGA-PBE 
eigenstates, whereas Go Wo on top of mPBEO overesti- 
mates the QP gap. A fully self-consistent GW calculation 
also yields a too large QP gap irrespective of the initial 
eigenstates, as a result of the neglect of the attractive 
electron-hole interaction^ By updating the eigenvalues 
(four times) in the Green's function G and keeping the 
screened Coulomb interaction W at the RPA level within 
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the initial PBE eigenvalues, the GWo@PBE scheme pro- 
duces a QP gap of 8.43 eV, in agreement with experi- 
ment. The VBM is now lowered by 2.77 eV with respect 
to the PBE eigenvalue, and the CBM is lifted up by 0.66 
eV, reaching to 4.90 eV by QP corrections. Compared to 
the GW result, the mPBEO CBM is placed 0.8 eV too 
high in energy, while the CBM energy calculated with 
the HSE06 functional coincides with that of the GWo 
(see Table |TlT|. A good agreement with the experiment 
F band absorption energy can be already obtained if the 
£ cbm hi Eq. (|2U|) is naively replaced by the GW value 
while keeping the ionization energy untouched. There- 
fore, the ionization energies of the neutral CI vacancy V° 
calculated by semi-local and hybrid functionals are well 
described from the total energy difference method. In 
contrast, the calculated ionization energy of the negative 
charge system V~ is less satisfactory with GGA-PBE and 
is not much improved with the hybrid functionals based 
on the experimental F' band absorption peak and the 
GWo CBM energy. This is easily understood since the 
electronic correlations for the removal of a second elec- 
tron from the a\ g level is beyond the scope of DFT. 91 
These many-body effects are accessible from the many- 
body perturbation theory, e.g. in the GW approximation 
via the self-energy. 

In the GW approximation, the QP energies of the high- 
est occupied and lowest unoccupied level correspond to 
the electron removal and addition energies, respectively. 
It is then straightforward from Eq. (|2U)) that the absorp- 
tion energy can be calculated as the QP energy differ- 
ence between the CBM and the lowest unoccupied state 
of the V + , provided the ionization potential of the neu- 
tral defect system is equivalent to the electron affinity of 
the positive charged system. For example, using a 64- 
atom supercell with a cutoff energy of 100 eV and a V 
point for the response function and 1024 total bands, the 
GWo@PBE method yields an F band absorption energy 
of 2.47 eV. We note that the two-particle excitonic effect 
in the optical absorption is not taken into account in the 
GW approximation either. 

For shallow defects, it has been found that the hybrid 
functional shows great improvement over local or semi- 
local functionals in the defect transition energies^ 8 . This 
is mostly likely benefited from the fact that the position 
of the shallow defect follows the band edge (either CBM 
or VBM), which can be reproduced by hybrid functionals 
with tunable a. For deep levels as demonstrated in this 
study, we find that a reproduction of a realistic band gap 
by an ad hoc tuning of the amount of the exact exchange 
in hybrid functionals does not guarantee an accurate de- 
scription of the optical defect levels, and the thermody- 
namic charge transition levels as well. We see that the 
GGA-PBE is prone to an underestimation of the vertical 
transition energy, which is obviously impaired by a small 
band gap and a low conduction band edge. A significant 
shift for the band edges can be observed with the hybrid 
functionals. This leads to an increased vertical transition 
energy which is usually in better agreement with exper- 



iment, although the overestimation of the CBM energy 
by an increased fraction of the exact exchange in mPBEO 
might give too high values (e.g. for the F band). 

In addition to the electronic contributions (e.g. ion- 
ization energy and the position of the band edge), the 
structural relaxation also plays an important role in pre- 
dicting the transition energies. The Frank-Condon shift 
(i.e. the difference between the absorption and emission 
energy) sheds light on the effect of the exact exchange on 
the lattice relaxations around the vacancy. By adopting 
the GWo CBM energy into the F band emission energy 
in Table [V] we find that the lattice relaxation in pres- 
ence of the electron-phonon interaction is best accounted 
by the mPBEO hybrid functional as a result of the more 
localized electron density in the vacancy. The localiza- 
tion is proportional to the amount of the non-local ex- 
act exchange which reduces the self-interaction arising 
from the DFT XC functional. The localized nature of 
the trapped electron is further enhanced by a GW cal- 
culation, exhibiting an even smaller k-dispersion of the 
dig level than the mPBEO since the GW approximation 
is free of self-interaction. In this context, we expect the 
structural relaxation will also be more realistic within 
many-body perturbation theory. 

The discussion of the vertical transition energies finally 
invokes us to return to the thermodynamic transition 
level presented in Sec. IIV CI The thermodynamic transi- 
tion level can be readily decomposed into the electronic 
and structural contributions as 

e(+/0) = [E g - E a (F)] + AE+ (21) 

and 

e(0/-) = [E g E a (F')} + A£° (22) 

where AE^ and AE^, are the relaxation energies from 
the initial geometry for the positive and neutral charge 
state, respectively. By incorporating the experimental 
vertical absorption energy and mPBEO relaxation energy 
into Eq. (|2"Tj) and (f2"2"f . we come to the thermodynamic 
transition energies e(+/0)=6.6 eV and e(0/— )=7.5 eV. 
Therefore, the general picture of the CI vacancy energet- 
ics by the adjusted hybrid functional so as to reproduce 
the experimental band gap remains qualitatively sound 
albeit not numerically accurate. 



VII. SUMMARY 

In this work we revisit the color centers in NaCl using 
hybrid density functionals. The reduced self-interaction 
error alleviated by the exact exchange in hybrid function- 
als allows for a more accurate description of the atomic 
relaxations in the vicinity of the vacancy compared to 
the DFT results. Yet, hybrid functionals are unable to 
achieve quantitative agreement with experimental opti- 
cal absorption peaks. We show that this is closely related 
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to the overestimations of the band edges when the funda- 
mental gap is reproduced by an empirical amount of ex- 
act exchange. Meanwhile, when the electronic correlation 
comes into play during the removal (addition) of a second 
electron from (to) the localized defect level, the ionization 
(affinity) energy predicted within the framework of DFT 
is far from satisfactory. More elaborate methods (such as 
GW approximation^-^ and two-particle Bethe-Salpeter 
equation^— ) are thus necessary for more accurate de- 
scriptions of the optical process and defect energetics of 



the localized defects in wide gap insulators. 
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